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C , Abstract. The Bose glass (BG) phase is the Griffiths region of the disordered Bose 

C/3 ' Hubbard model (BHM), characterized by finite, quasi-superfiuid clusters within a Mott 

insulating background. We propose to utilize this characterization to identify the 
complete zero-temperature phase diagram of the disordered BHM in rf > 2 dimensions 
by analyzing the geometric properties of what we call superfiuid (SF) clusters, which 
are defined to be clusters of sites with non- integer expectation values for the local boson 
'^ I occupation number. The Mott insulator (MI) phase then is the region in the phase 

fl ' diagram where no SF clusters exist, and the SF phase the region, where SF clusters 

percolate - the BG phase is in between: SF clusters exist, but do not percolate. This 

definition is particularly useful in the context of local mean field (LMF, or Gutzwiller- 

Ansatz) calculations, where we show that an identification of the phases on the basis 

^ ] of global quantities like the averaged SF order parameter and the compressibility are 

0^ ' misleading. We apply the SF cluster analysis to the LMF ground states of the two 

L^ . dimensional disordered BHM to produce its phase diagram and find a) an excellent 

If-y I agreement with the phase diagram predicted on the basis of quantum Monte Carlo 

simulations for the commensurate density n ~ I, and b) large differences to stochastic 
f~^ , mean field and other mean field predictions for fixed disorder strength. The relation 

CO ' of the percolation transition of the SF clusters with the onset of non-vanishing SF 

stiffness indicating the BG to SF transition is discussed. 
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1. Introduction 

The experimental proof of tlie Mott insulator (MI) to superfluid (SF) transition in 
ultracold atomic systems [1] opened a wide field of interesting research in this field. In 
particular the influence of disorder on a system of bosons in a regular (e.g. optical) 
lattice received much interest since the fundamental work of Fisher et al. [2] . Here the 
phase diagram and transitions for bosons in a disordered potential was analysed and 
the existence of a Bose glass (BG) phase was predicted. The BG represents a non- 
SF but, in contrast to the MI, compressible phase displaying an excitation spectrum 
with arbitrarily small excitation energies. The BG phase is the analogue of the Griffith 
regions occurring for instance in disordered magnets, whose physics is dominated by 
rare region effects due to arbitrarily large strongly coupled clusters [3, 4, 5]. 

Studies of the excitation spectrum of the disordered system in dependence of the 
disorder strength and time-of-flight measurements confirmed the predicted BG phase 
experimentally [6]. In addition to the well controllable optical lattice, disorder is 
introduced either by a non-commensurate periodic potential [6] or by speckle potentials 
[7, 8]. A new view on the properties of ultra cold bosonic gases opened up as high 
resolution techniques allowed access to single-site detection recently [9, 10]. This 
progress now yields a direct view on the population numbers within the different phases 
and the in-situ hopping dynamics of the bosons in their optical potential. 

Theoretically the phase diagram of the disordered Bose Hubbard model (BHM) has 
been studied by various methods: The strong coupling expansion [11] is a perturbative 
method up to third order in the tunnelling rate yielding a prediction on the Mott lobes. 
The disordered BHM was widely studied by quantum Monte Carlo (QMC) methods in 
various dimension [12, 13, 14, 15, 4, 16, 17]. In addition, density-matrix renormalization 
group (DMRG) techniques were applied to ID systems containing either quasi-periodic 
potentials [18, 19, 20] or a uniform distribution of disorder strength [21, 22]. A frequently 
used alternative approach is the local mean field (LMF) approximation [23], which 
replaces the nearest neighbour hopping on the lattice by isolated bosonic degrees of 
freedom interacting via an effective mean field coupling with the neighbours. Based 
on the LMF approximation several numerical techniques, such as stochastic mean field 
(SMF) theory [24, 25] and LMF theory [26, 27, 28], were proposed. 

An intriguing question has been for a long time the potential existence of a direct 
MI-SF transition [12, 13, 14, 15, 4, 16, 17], which is now excluded by a rigorous theorem 
[13]. The occurrence of the BG phase intervening between the MI and SF is caused by 
Griffiths effects [29, 2] due to arbitrarily large, but exponentially rare clusters of one 
phase within a background of another phase [30, 31, 3]. Since any exponentially rare 
event is hard to sample numerically, the existence of an intervening BG phase might 
have been eluded some studies, be it QMC [14, 15, 4, 16, 17], LMF [26, 27, 28] or SMF 
theory [25, 24]. One might speculate that, rather than calculating spatially averaged 
quantities, a look at the aforementioned clusters in individual disorder realization itself 
would tell us more about the actual state the system is in. In this paper we propose a 
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method to identify the different phases of the disordered BHM on the basis of geometric 
properties of what we call SF clusters, which are clusters of sites with non-integer boson 
occupation number. The MI phase then is the region in the phase diagram where no SF 
clusters exist, and the SF phase the region, where SF clusters percolate - the BG phase is 
in between: SF clusters exist, but do not percolate. We apply this criterion to results of 
LMF calculations and compare it with predictions of other methods: on one side SMF 
theory, where the individual phases are identified on the basis of spatially averaged 
LMF quantities like SF order parameter or compressibility, and on the other side QMC 
simulations, which are supposed to be exact up to statistical and extrapolation {L — )■ oo, 
T — )■ 0) errors (which can, of course, be quite large). Our aim is to demonstrate that 
the use of averaged quantities in LMF theory leads to incorrect predictions and that the 
cluster analysis predicts the phase diagram in rf > 2 in very good agreement with the 
exact QMC results, even when applied to LMF data. 

The paper is organized as follows: In section 2 we recapitulate the LMF and SMF 
theory for the disordered BHM. In section 3 we critically examine the use of the averaged 
SF order parameter ip and the compressibility k as indicators of the different phases of 
the disordered BHM in LMF and SMF theory and then introduce our new method 
to construct the phase diagram on the basis of an analysis of SF clusters. We apply 
this cluster analysis in section 4 to the 2d disordered BHM with commensurate filling 
and with fixed disorder strength and compare it with prediction from QMC simulations 
and from SMF theory. Section 5 concludes with a discussion of the equivalence of the 
predictions of the SF cluster analysis for the phase boundaries with the conventional 
definition of the MI-BG and BG-SF transition points. 

2. The model 

Ultracold bosonic atoms in an two dimensional square optical lattice can be described 
by the BHM 



^ = J2^^i~^^')^i + ^J2^^\^i~^) "^Zl^^^^i' (^) 



where i = 1, . . . , M is the site index, M = L'^ the number of sites and L the lateral size of 
the square lattice. The chemical potential is described by fj,, the inter particle repulsion 
by U and the tunnelling rate by J. The last sum runs over all nearest neighbour pairs 
(ij) of the underlying lattice. The operator n^ = ala^ is the particle number operator of 
bosons on site i, which are annihilated and created by the operators a- and a|. Moreover, 
on-site disorder is introduced by the parameter e,, which is drawn randomly from a box 
distribution p{e) = O (A/2 — |e|) /A, where A is the strength of the disorder. 

In the SF regime tunnelling dominates the system. Thus, the ground state is a 
coherent state, which is an eigenstate of the tunnelling part of the Hamiltonian. The 
SF parameter ipi = (dj) is the expectation value of d^ evaluated in the ground state for 
T = 0, which is non-zero for a coherent state. An eigenstate of the diagonal part of the 
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Hamiltonian on the other hand is a Fock-state, which is the ground state of the system 
for small tunnelling rates in the MI regime. The expectation value of a^ in a Fock-state 
is zero in any case. Hence, the mean value of the SF parameters ifj = ^^ ipi/M is an 
order parameter for the SF phase. 

The phase diagram of the pure system in dependence of JZ/U and /x/f/ consists 
of so called Mott lobes, in which the system is a MI with a fixed integer number of 
atoms per site. While in the MI regime the state is localized in real space, in the 
surrounding SF regime it is localized in fc-space. When disorder is introduced, the Mott 
lobes shrink and a new phase, the BG phase, occurs. In order to distinguish all three 
phases the compressibility k = {n^) — (n)^ is necessary. Among the three phases only 
the MI is non-compressible and only in the SF phase the SF order parameter is non- 
zero. Consequently, the phase of the system can be identified by the SF order parameter 
ip = Y2^ ipi/M and the compressibility k. While the SF order parameter is a measure for 
the coherence in the system, the compressibility describes the variance of the particle 
number per site. 

In order to determine the phase diagram of the BHM, the ground state properties 
of the LMF Hamiltonian can be studied via SMF theory [24, 25], which computes the 
PD self-consistently, or via LMF theory [23, 28, 27], which solves the coupled set of 
equations for the local SF parameter directly on the lattice. We will first describe 
the approximations made in LMF theory, followed by a discussion of the additional 
assumptions made in the SMF approach. 

2.1. Local mean field theory 

In LMF theory the tunnelling part of the Hamiltonian can be approximated via 

a^a] ^ a^{a]) + a^a^) - (a,) (a]), (2) 

where terms of the form (a^ — {a^)){a}j — {oj)) are neglected [23]. The central quantities 
are the local SF order parameters 

V'i = {9s\ai\gs), (3) 

which are defined as the expectation values of the annihilation operator at the individual 
site i in the ground state of the system. Because of the [/(l)-symmetry they can be 
chosen to be positive and real, which leads to 

a^a^j ~ i^jO-i + V'jaj - i'ii'j. (4) 

Thus, the Hamiltonian can be decomposed in a sum of diagonal operators, 

i 

Hi = {e^ -fi)% + —fii (n. - l) - Jr}i («i + «1 - ^i) , (5) 
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whose tunnelling rate is replaced by an effective local rate Jr]i, which depends on the 
local SF parameter of the neighbouring sites rji := Ylj^ij'^j^ with Aij = 1 for i and 
j nearest neighbours on the square lattice with periodic boundary conditions and zero 
otherwise. This approximation reduces the full quantum problem to M quantum sites, 
which are coupled in a mean field way with a spatially varying coupling rate. 

In order to compute the phase diagram in LMF theory the coupled set of the self- 
consistency equations 

^, = {a,), z = l,...,M, M = L^ (6) 

is solved on a LxL lattice, where the expectation value is evaluated in the ground state of 
H^ that itself depends an the local SF parameter ipi. As a result of the decomposition (5) 
all states considered (in particular the ground state) are a direct product of individual 
single-site states. This means in particular that they are Gutzwiller states of the form 

M / oo \ 

i=l \n=0 / 

with single-site states |0j) = J2'^=o^n\^)i given in particle number basis and |c^p 
describes the probability to find n bosons at site i and fulfils ^^^^ |c^P = 1. In 
particular, the local SF order parameter is then given by ipi = ^^o^n-i*^™^/^ ^^d 
the local boson number is represented by (h-) = ^^q I'^nP^*- 

For the numerical implementation, starting from a random initial configuration 
for tpi on the 2d lattice, the set (6) of equations is solved recursively. This involves 
solving the eigenvalue problem on each site and computing the expectation value of 
the annihilator in the numerically determined ground state. This is repeated until the 
averaged SF order parameter 



V' 



1 ^ 
M ^^ 



(8) 



is determined with an accuracy of 10~^. In the disordered case the results are averaged 
over 200 different realizations of disorder, indicated by the brackets [. . .]av Since we are 
working in a regime in which the maximum average particle number per site is three, it 
was numerically checked that it is sufficient to truncate the basis of the Hilbert space 
for each site at iV = 10. With the solution found for the local SF parameters ipi on the 
lattice, the ground state of Hamiltonian (5) is calculated numerically. Afterwards all 
desired expectation values and finally the compressibility 

(iV2) - (iV)2 

with A^ = Y2i f^i can directly be computed. The probability distribution (PD) 

1 ^'' 



(9) 



(10) 



is determined on the basis of the complete set of values of ipi, additionally averaged over 
disorder realizations. 
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2.2. Stochastic mean field theory 

The central idea of SMF theory is to circumvent the computation of all local order 
parameters ipi by deriving a self-consistency equation for the probability distribution 
P{ip) directly [24, 25]. Additional approximations are of course necessary. The SF order 
parameter ijj = {gs\d\gs), which is derived from the ground state of the full quantum 
Hamiltonian, can be determined by the ground state of the single-site Hamiltonian (5) 
in dependence of the stochastic variables e and rj. The parameter e is then a stochastic 
variable drawn from the disorder distribution p{e) and as a result ip = {gs\a\gs) 
is a stochastic variable, drawn from the PD P{ip), which must be determined self- 
consistently. Since r] is the sum of the SF parameters of the neighbouring sites it also 
is a stochastic variable drawn from Q (77). The problem of computing the ground state 
of the full quantum system for all lattice sites simultaneously is thereby replaced by 
analysing the ground state \gs (e,?])) of the site independent Hamiltonian 

H = {e- n)h-\ n {n - 1) ~ Jt] {a + a^ - tp) (11) 

as a function of e and 77. Thus, the probability distribution (PD) 

Pi^) = jdriQiv)Pvi^) (12) 

depends on the distribution Q (rj) of the occurring values of rj and the distribution P^ (ip) 
of the local SF parameters for given rj. A direct analysis of {gs{e,ri) \d\gs {e,ri)) as a 
function of e and f] yields 

Since rj is the sum of the local SF parameters ip of the neighbouring sites its distribution 
is given by 






(14) 



where P^ (^/'i, . . . , ^Z'^) is the connected probability distribution function of the local 
order parameters t/^i, . . . ,^z of the Z neighbour sites of a single-site. Assuming that 
these Z local SF parameters are statistically independent 

z 

Pz(v^i,...,^z)=n^(V'0, (15) 

equation (14) transforms into a convolution 



,i=l / \ i=l 
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Since the assumption (15) implies the absence of correlations of the local SF parameters, 
one expects that it is not justified close to the phase boundaries, where the correlation 
length even diverges, when the transition is 2nd order. We examine the validity of this 
approximation in dependence of the system parameter JZ/U and ^/U in Appendix A. 
After determining the PD P [ip) the SF order parameter ip = J dtpP (tp) ip is given 
by the mean value of the distribution. The compressibility k = [{N"^) — (iV)^]av with 
N = Y2i ^i is computed. With these quantities at hand one can, on the basis of the 
underlying approximations, estimate the phase boundaries of the transitions between 
MI, BG and SF. They are shown in Figure 4 and 5 and will be discussed in the following. 

3. Criterion for the phase transition 

In this section we will first discuss the well known SF order parameter ip and the 
compressibility k, which are expected to indicate the phase transition: The ground 
state in the MI regime is a Fock-state, which is incompressible (k = 0) and non-coherent 
z/^ = 0, while conversely in the SF regime it is described by a coherent state ■?/' = 0, which 
is compressible (k > 0). If disorder is introduced, those phases are separated by the BG 
phase, which is compressible (k > 0), but not coherent V' = 0. We will see that a precise 
prediction of the transition point on the basis of ^ or k is not possible in LMF theory. 
Instead we will introduce an identification criterion of the different phases on the basis 
of the complete set of local occupation numbers. 

3.1. SF order parameter and compressibility 

In the ordered case the on-site energies e^ are zero and the lattice is homogeneous. The 
SF order parameter ip clearly marks the location of the MI to SF phase transition as 
shown in Figure 1 on the left hand side for ix/U = 1.05 and 0.32. While the SF order 
parameter is zero for small tunnelling rates in the MI phase, it becomes non-zero and 
positive above a critical value of the tunnelling rate in the SF phase. The compressibility 
K shows the same behaviour at the phase transition as the SF order parameter for both 
methods in the ordered case. Moreover, we analysed different lattice sizes L and the 
LMF results show no visible finite-size effects. In this way the phase transition in the 
ordered case can be determined very precisely, both within SMF and LMF theory. The 
resulting phase transitions agree perfectly with the perturbation predictions [32] 

^ = -^f -2- + l)±t/7(^-l) -"-^n, (17) 

where n denotes the mean number of particles per site and simultaneously counts the 
number of lobes. The calculation in [11] predicts that in the disordered case the upper 
(lower) part of the Mott lobes are shifted downwards (upwards) by A/2 but the shape 
remains unchanged. 
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Figure 1: Comparison of the LMF and SMF predictions for the average SF order parameter -0 and 
the compressibility k for fixed chemical potential /i and varying tunnelling rate J. Left: Homogeneous 
case (A = 0), Right: disordered case with A/U ~ 0.6. Top row is for jj, = 1.05, where the ordered 
system displays a MI SF transition and the disordered system a BG SF transition (k > for all values 
of J). Bottom row is for fj, = 0.32, where the ordered system again displays a MI SF transition and 
the disordered system is expected to display MI, BG and SF phases (see section 3.2). For LMF theory 
the results for a 2d lattice with L = 100 (line), L = 50 (o), L = 10 (+) are depicted, which shows that 
finite-size effects can be neglected. 



The situation for the disordered case is shown in Figure 1 on the right hand side. 
The SF order parameter is shown for fi/U = 1.05 and 0.32 as a function of the tunnelhng 
rate for a disorder of A/U = 0.6. Whereas the SMF theory predicts a direct BG SF 
transition, at a critical value JZ/U ~ 0.0241 and 1.0455, see 1 (b, d), above which ijj 
become non-zero the behaviour of ip as predicted by LMF theory does not indicate a 
transition at all; it varies smoothly with the tunnelling rate J. This is not a finite- 
size effect as we have checked by examining different lattice sites, as shown in 1. The 
compressibility, which indicates the MI BG transition, displays the same behaviour. 

It turns out that the reason for the failure of the average SF order parameter 
to predict the location of the BG SF boundary is the following: In the disordered 
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case the value of the local SF parameter varies substantially from site to site due 
to the variation of the local potential of e^. Close to the phase transition there are 
sites with zero local SF parameter and others, where the local SF parameter is still 
positive. This has been interpreted as an overestimation of the phase coherence in LMF 
description [24]. Our interpretation, however, is different: It is only the average SF order 
parameter ip that overestimates the phase coherence. A closer look at the complete 
probability distribution P(^) of the local SF order parameters and their geometrical 
features provides an estimate of the SF regions in the phase diagram. Its prospects are 
discussed in the 4.3. In the next section we discuss, how a deeper understanding of the 
mechanisms driving the phase transitions and their location in the phase diagram can 
be obtained by studying the geometric characteristics of the spatial inhomogeneities of 
the local SF parameters ipi and particle number per site (n^). 

3.2. Identification of phases via local boson occupation number 

The MI and SF phases can be discriminated by the boson number statistics at individual 
sites, as has also been demonstrated experimentally in [1]. The ground state in the 
extreme MI limit (J — )■ 0) is a Fock states with a definite number of particles n at 
each site. In the extreme SF limit {U -^ 0), the ground state is a coherent state , in 
which the local boson number distribution is close to a Poissonian. Although, in the 
regime between these two extreme limits the ground state wave function can no longer 
be written as simple product states still the MI phase is characterized by a sharp, integer 
boson number per site and the SF phase by a fluctuating boson number per site, i.e. a 
no n- vanishing variance of the boson number distribution pl^ = |c^p (c.f. the expansion 
coefficient in the Gutzwiller wave function (7)). In other words in the MI regime the 
expectation value of the number of bosons per site (n-) is an integer, whereas in the SF 
regime it is non-integer. 

Whereas in the ground state of the homogeneous BHM either all sites have an 
integer boson number (MI regime) or all sites have a non-integer boson number (SF 
regime) this is different in the disordered BHM. In particular, outside of the MI regime 
one expects to encounter spatially inhomogeneous situations, in which some sites have 
a sharp (integer) boson occupation numbers and others have fluctuating (non-integer) 
boson numbers. Introducing phase operators $ that is canonically conjugate to the 
boson number operators fii the BHM can be mapped upon a Josephson junction array 
or more general to a quantum rotor model [2], in which superfluidity is indicated by long 
range order in these phases {d >2). Because of the Heisenberg uncertainty relation sites 
with sharp phases correspond to sites with fluctuating boson numbers, and connected 
clusters of sites with fluctuating boson numbers tend to have, roughly speaking, all the 
same phase. These clusters can therefore be identifled with SF regions, although, true 
superfluidity only exists in the inflnite system. Indeed, once these phase ordered clusters 
percolate, true superfluidity emerges, signifled by a non-vanishing SF stiffness, which is 
the extra free energy cost to impose a uniform twist on the phases. Since such a uniform 
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twist can be introduced by applying a certain twist at the boundary phases in one space 
direction, it is clear that the SF stiffness is zero as long as the clusters do not percolate: 
In the absence of long range order in the phases such a twist at the boundary does not 
cost energy. 

On the basis of this qualitative picture we hypothesize that the BG to SF transition 
in the d-dimensional BHM [d > 2) coincides with the percolation transition of connected 
clusters of lattice sites with a non-integer boson number expectation value (nj. We 
expect this coincidence to hold as long as the SF phase displays true long range order, 
characterized within the phase description by a non-vanishing long distance limit of 
phase correlations - which means it should hold for d > 2. The BG SF transition in one 
dimensional BHM might not be related to a percolation transition, since in d = 1 the SF 
phase has only quasi long range order (algebraically decaying correlations). We should 
note that the relation between the BG SF transition and percolation has already been 
pointed in [33, 26, 34], but has neither been used in a quantitative manner to determine 
phase boundaries nor checked against, for instance, Monte Carlo results. 

In the following we denote the sites with non-integer boson occupation number (fi^) 
as SF sites, and sites with integer (nj as MI sites. Analogously, we discriminate SF 
clusters and MI clusters. Formally we map the boson occupation numbers to a discrete 
field Si that is set to Si = 1 for SF sites and Si = for MI sites. Then we identify the 
different phases of the disordered BHM as follows: 

MI phase: Si = for all sites i. All boson occupation numbers are integer (and 
identical), consequently the compressibility k, is zero. 

SF phase: Si = 1 for a macroscopic fraction of sites, which form a percolating 
connected cluster. According to what we discussed above the percolating cluster has 
phase long range order and thus yields a non- vanishing SF stiffness (which is proportional 
to the SF density). 

BG phase: Characterized by a non- vanishing density of sites with Si = 1, none 
of the connected clusters formed by the SF sites percolates. The BG phase is thus 
characterized by isolated SF clusters within a MI sea. The phases of the isolated clusters 
are uncorrelated, hence phase long range order is lacking and the SF density vanishes 
(no SF order). Moreover, due to the number fluctuations on the SF sites the BG phase 
is compressible (k > 0). 

Within LMF theory the expectation values of the local boson occupation numbers 
are straightforward to calculate via rij = {gs\n^\gs), wheTe\gs) is the ground state of 
the LMF Hamiltonian (5). For numerical reasons we introduce a threshold 7„ into the 
definition of the discrete field 

. ,- - - ln<{n,)<I + ln / = 0,1,2, 




where 7n = 5 ■ 10"'^ is chosen to serve as the cut-off in this algorithm. In the whole 
parameter range, where sites with integer particle number occur, the histogram of 
the mean particle number (rij) has narrow peaks of width 7„ at integer values. The 
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Figure 2: Configurations of tlic local SF parameter ipi, the occupation number {fi^} and the discrete 
variable Si (18) for a single realization of disorder for A/U = 0.6. The first row shows an example for 
MI {JZ/U = 0.0242, n/U = 0.4394) followed by one for BG {JZ/U = 0.0182, ^/U = 1.0455) and SF 
{JZ/U = 0.141, ^/U = 1.0455). Note that blue marks the minimal value (zero in the left and right, 
one in the middle column) and green the maximal value (one in the left and right, two in the middle 
column) . 



width decreases when we increase the number of iteration steps to solve the self- 
consistency equations 6. The threshold parameter 7„ introduced to identify MI sites (and 
complementarity SF sites) can be reduced by increasing the numerical effort without 
changing the final results. 

In Figure 2 typical results for one realization of disorder for A/U = 0.6 and 
fi/U = 1.0455 are shown for three different values of the tunnelling rate JZ/U for 
the three phases . In the first row the local SF parameter ipi, in the second the particle 
number per site (n/) and in the third the resulting discrete map Si is shown. In the MI 
regime all sites are occupied by the same integer number of particles (in this case one, 
since we are in the first Mott lobe). At the transition from the MI to the BG regime SF 
sites {Si = 1) with non-integer particle number occur. Because of these locally occurring 
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SF sites ipi > the SF order parameter ip is small but not zero in this regime. Since 
the SF islands are compressible, this phase has positive compressibility. In the BG 
phase the SF islands does not percolate, yet. They grow in number and size, until 
one of them finally percolates. The percolation represents the actual transition to the 
SF regime in parameter space. Just after the percolation the phase in the system is 
coherence macroscopically, which means that all local SF parameter ipi are positive and 
compressible, as described above. 

3.3. Percolation analysis 
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Figure 3: The percolation probability ppcic of the SF cluster (top) and finite size scaling plot 
(bottom). The critical tunnelling rate according to the finite-size scaling (19) are given by JcZ/U = 
0.15 for ^i/U = 0.439 (left) and JcZ/U = 0.04 for ^i/U = 1 (right); the critical exponent is i^ = 1.33. 



In this section we demonstrate how we determine numerically the percolation 
transition via a cluster analysis of the discrete map Si and finite size scaling [35]. 
Assume we study the phase diagram in dependence of the system parameter called 
X and y. Than, the percolation probability ppcrc? i- e. the probability of having a 
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percolating cluster, is given for fixed y as a function of x and will be determined for 
different system sizes L. The percolation probability is expected to obey the finite-size 
scaling form 

pp,,,{L,x)=p{L'/''{x-x,)), (19) 

where Xc is the percolation threshold, i.e. the value above which a percolating cluster 
exists with probability one, and u the critical exponent determining the divergence of 
the mean lateral cluster size at the transition. The scaling function p{X) approaches 
zero for X <^ 1 and one for X ^ 1, which means that exactly at the transition Xc the 
curves for different system size should intersect (in the scaling limit). This intersection 
point, which we can easily be identified with the system sizes behaviour at hand, is thus 
a reliable indicator for the percolation transition. 

The cluster analysis of the discrete map Si is done for every disorder realization. 
Afterwards the results are averaged over 200 {L = 50, 100) and 2500 [L = 10) 
realizations of disorder. The percolation probability ppcrc for this case is shown in 
Figure 3 (blue for MI and red for SF sites) for different system sizes as a function of 
the tunnelling rate JZ/U. Moreover, the finite size scaling analysis for the percolation 
transition at fi/U = 0.439 and ^/U = 1 for A/U = 0.6, yielding JcZ/U = 0.15, 
JcZ/U = 0.04 respectively and the critical exponent i^ = 1.33 in both cases, is depicted. 
Thus, this transition is in the universality class of conventional 2d percolation [35]. We 
find the same universality class of the percolation transition for all parameter values 
that we studied. 

4. Results 

4-1. Commensurate filling - Comparison with QMC results 

In this section we determine the complete phase diagram for commensurate density in 
dependence of A/2 J and U/J, for which a prediction on the basis of QMC simulations 
is available [12]. We fix the particle density to n = (^^=1^^) / M = 1 with an accuracy 
of 10~^ by adjusting the chemical potential for each point (A/2J, U/J) in the phase 
diagram that we study. Outside of the Mott lobes this result is unique, whereas in the 
MI regime the chemical potential is fixed to the middle of the MI gap. In the fi/A 
versus JZ/U representation, where the Mott lobes are visible and which we will discuss 
in section 4.2, this n = 1 line always passes the tip of the first Mott lobe. In the A/2 J 
versus U/J parameter space the corresponding line for fixed A is a straight line through 
the origin with slope A/2U. 

With the chemical potential that fixes the density n to one we compute the ground 
state of the LMF Hamiltonian (5) and determine the discretized boson number field Si 
(18), which we use to identify MI, BG and SF phase. The resulting phase diagram is 
shown in Figure 4 on the left. As expected [13] the SF region is completely surrounded 
by the BG phase (except at A = 0). Its boundary has some characteristic features: It 
extends in a slight bump up to quite large disorder strength (up to A/2J ~ 75) and 
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in a pronounced nose up to the interaction strength U/J ~ 52. This nose gives rise 
to a re-entrant behaviour: Moving vertically from a point within the MI phase, which 
has long range positional order, one enters first the BG phase, which is disordered and 
then, upon further increasing the disorder strength, enters the SF phase, which has 
off-diagonal long range order. Weak disorder thus supports superfluidity in the BHM, 
as has been observed before [17, 12, 36, 37]. 
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Figure 4: Left: LMF cluster analysis phase diagram for commensurate boson density n = 
(^j^-^nj)/A/ = 1 determined with the discretized boson occupation number field Si (18). The 
percolation transition of the SF sites (Si = 1) occurs when crossing the red line, which indicates 
the BF-SF phase boundary. The blue line marks the boundary of the MI region, in which all sites are 
MI sites (Si = 0). The black line indicates the MI-BG transition according to the perturbative result 
(17). Right: Prediction for the phase diagram for commensurate boson density n — 1 based on the 
results of QMC simulations (data taken from [12]). 



Remarkably, our prediction on the basis of a cluster analysis of LMF ground states 
agrees very well with the results of QMC simulations [12] shown for comparison in Figure 
4 on the right. The shape of the SF-BG phase boundary with its characteristic nose 
and bumps clearly coincide. The quantitative agreement is very good, too, regarding 
the substantial error bars of the QMC data in the large disorder and large interaction 
regime (the QMC estimate for the extreme value of A in the bump is A/2 J ~ 72 ± 4 
and of inter particle interaction in the nose U/J = 49 ± 3, c.f Fig. 2 in [12]). Moreover, 
with our method we could also explore the weak interaction region, which is hardly 
accessible by QMC methods, and found a singular behaviour of A with [/ — )■ 0, which 
is compatible with the analytically predicted behaviour A oc \/U [38]. We conclude 
that the percolation criterion that we introduced in section 3.2 to locate the SF-BG 
transition produces remarkably accurate predictions even in LMF theory. 

Our result for the MI-BG transition line, which denotes the appearance of non- 
integer boson occupation numbers and thus SF sites, agrees well with the perturbative 
result 17, shown in Figure 4 on the left. Moreover, they agree with the line A = Eg/2 
obtained using the gap data from [39], shown in Figure 4 on the right. 

In passing, we note that for weak disorder the MI clusters percolate close to the 
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BG-SF transition line, whereas for stronger disorder they percolate deeper inside the 
BG phase. Whereas for weak disorder the individual sites of a MI cluster in the BG 
phase all have the same integer occupation number, this is in general not the case any 
more for strong disorder: the integer occupation number of MI clusters can vary from 
site to site. 

Finally, we note that SMF theory as described in section 3.1 predicts a direct MI- 
SF transition along the lower border of the SF region in the parameter range shown 
in Figure 4. The characteristic BG region for small disorder strength is absent in this 
parameter range, which is in contradiction to the theorems proven in [13], which exclude 
a direct MI-SF transition in any disordered systems. Besides, we checked that the BG 
phase occurs for even higher values of U/J, which is in agreement with results to be 
presented in the next section. 

4-2. Fixed Disorder - Comparison with SMF theory 

After we have seen in the last section that our method to determine the phase diagram of 
the 2d disordered BHM leads to results that agree very well with QMC predictions, we 
determine in this section the [fi/U-JZ/U) phase diagram for a fixed disorder strength 
A/ J = 0.6 and compare it with predictions of SMF theory. In this phase diagram the 
Mott lobes occur and the hne given by (n) = 1 always passes the tip of the first one. 
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Figure 5: Comparison of LMF cluster analysis and SMF phase diagram for fixed disorder strength 
A/U ~ 0.6. Left: LMF cluster analysis phase diagram determined with the discretized boson 
occupation number field Si (18). The percolation transition of the SF sites {Si = 1) occurs when 
crossing the red line, which indicates the BG-SF phase boundary. The blue line marks the boundary of 
the MI region, in which all sites are MI sites {Si = 0). The black line is the MI-BG transition according 
to the perturbative result given by (17). Right: SMF phase diagram determined by using the SF order 
parameter tp and the compressibility k [24, 25]. The red line indicates the critical tunnelling rate J 
where the SF order parameter ip becomes non-zero, the blue line the critical tunnelling rate, where the 
compressibility k becomes non-zero. 



In section 2.2 we introduced SMF theory and already emphasized that SMF theory 
bases on the same approximation to the Hamiltonian as LMF theory, but, it involves 
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MI-BG 

JZ/U fi/U 


BG-SF 

JZ/U fi/U 


A/U 


LMF (cluster analysis) 
QMC results [12] 
strong-coupling expansion [11] 


0.1115 0.4644 
0.124 (0.4561) 
0.1253 0.4345 


0.1509 0.4434 
0.2012 (0.4082) 


0.6 
0.6 
0.6 


LMF (cluster analysis) 
QMC results [12] 




0.0942 0.4868 
0.1047 (0.4846) 


1 
1 


LMF (cluster analysis) 
QMC results [12] 
quantum rotors model [37] 




0.0934 0.5043 
0.1062 (0.4950) 
0.112 0.375 


2 
2 
2 



Table 1: Comparison of the parameters at the tip of the first Mott lobe of quantmii models (1) with 
LMF cluster analysis results. The chemical potential /i/C/ for the QMC results of [12] is our LMF 
estimate for a density n = 1 and fixed values of JZ/U and A/2 J. The BG-SF predictions of [12] were 
not obtained by QMC of the disordered BHM, but are based on gap data for the ordered BHM [39]. 



the additional approximation (15) on the distribution Pz {tpi, . . . ,ipz)- The validity of 
this restriction fails close to the phase transitions as we show in Appendix A. Despite 
or perhaps because of this approximation the SF order parameter ip as well as the 
compressibility k, computed within SMF theory are exactly zero in specific regions of 
the parameter space (c.f. Figure 1), which one might want to identify with MI and EG 
phase, as was done in [24, 25], c.f. section 3.1. The LMF cluster analysis and SMF 
[24, 25] phase diagrams for fixed disorder strength A/U = 0.6 are shown in Figure 
5. One immediately observes substantial differences: Firstly, in LMF theory the BG 
phase always separates the MI from the SF phase. The intervening BG phase is actually 
predicted to be quite large even at the tip of the Mott lobes, not just a "thin sliver" 
[2]. SMF theory, however, predicts a direct MI-SF transition, in contradiction to [13]. 
Secondly, large differences in the critical tunnelling rate for the BG-SF transition occur 
especially in the region around /i = 1. Assume we fix the chemical potential there. 
In this case the SMF theory predicts the phase transition at JZ/U = 0.0241. The 
percolation of the SF cluster, however, takes place at JZ/U = 0.0430. Thus, significant 
changes of the system in this case occur for values of the tunnelling rate twice as large 
as predicted by ip in SMF theory. 

A direct comparison of our results with the QMC data of [12] is not possible here, 
since the latter are obtained for the canonical ensemble, where the chemical potential 
is absent. However, we can take our LMF estimate for the value of /x that fixes the 
particle density at n = 1 for fixed U/J and A/2 J to obtain an approximate comparison 
- see Table 1, where we also show the prediction of the strong coupling expansion [11] 
for the MI-BG transition for A = 0.6. One observes deviations of the QMC and strong 
coupling predictions from our LMF cluster analysis results at the tip of the Mott lobe, 
but a good agreement for stronger disorder, A/U = 2. At this disorder strength also a 
QMC prediction for the quantum rotor model exists [37], which differs by 25% from the 
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predictions for the BHM, ours and the one from QMC We also note that the tapered 
shape of the Mott lobe predicted by the strong coupling expansion [11] agrees well with 
our result of the LMF cluster analysis shown in Figure 5. 

In addition to our LMF cluster analysis and the SMF theory discussed above 
a number of other approximative methods have been applied to calculate the phase 
diagram of the BHM at fixed disorder. In the zero temperature mean field phase diagram 
of [40] the EG phase is completely absent, which might be true in infinite dimensions, 
but certainly not in c? = 2 or 3. 

A LMF theory has been used in [27, 26] to solve the self consistency equations (3) 
and to calculate a LMF expression for the stiffness or SF fraction and the compressibility. 
Using these two observables the (/x/f/- JZ/f/)-phase diagram is then determined, which 
displayed a round shape of the Mott lobes, a direct MI-SF transition for small disorder 
and an intervening BG phase at larger disorder. It should be noted that although the 
starting point of the calculation in [27, 26], the LMF approximation, is identical to ours, 
the usage of a different criterion to identify the phases leads to a phase diagram that 
differs significantly from the one predicted by us. 

A multi-site LMF theory is introduced in [41], where every plaque of two by two 
sites is treated quantum, which keeps the spatial correlation therein. Instead of single 
sites these plaques are coupled in a LMF way (analogous to section 2.1). The Mott lobe 
is determined for both, the single-site and multi-site LMF theory, on the basis of the 
condensate fraction. In agreement with [27] and SMF theory it shows a round shape 
at the tip. The multi-site LMF theory predicts a larger MI region than the single-site 
LMF theory. Note that the condensate fraction smoothly approaches zero, analogous 
to our observations on the SF order parameter and the compressibility made in section 
3.1; a linear fit is used to determine the transition point. 

In [42] the so-called the Gutzwiller projected variational techniques is introduced 
in order to determine a canonical transformation of the quantum Hamiltonian, which 
requires the truncation of the hopping term. Thus, it is possible to minimize the 
expectation value of the transformed Hamiltonian in Gutzwiller type local mean field 
states with respect to its variational parameters. Finally, the SF stiffness and the 
compressibility yield the phase diagram, which shows a remarkably narrow BG region 
between the MI and the SF phase. 

In all mean field calculations mentioned the tip of the Mott lobe is predicted for far 
higher values of the tunnelling rate as our results based on the LMF cluster analysis, 
and results of QMC or strong coupling expansion methods, as listed in Table 1. 

4.3. The probability distribution of the local SF parameter 

In this section we discuss the probability distribution (PD) of the local SF order 
parameter P{ip). In the ordered case A = 0, depicted in the first row, all local SF 
parameter are identical since all sites have the same on-site energy. The averaged order 
parameter ip, depicted as a blue cross, is identical to each local SF parameter ipi and 
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the variance of this value is zero. Therefore, the PD P (ip) is a delta function at the 
values of ip. Within the MI region the local SF parameter is zero everywhere and thus 
the PD is a delta function at ip = 0. In the SF regime still the PD is a delta function 
but at positive ip, which increases with the tunnelling strength. 
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Figure 6: Probability distribution P (-0) of the local SF parameter at fixed chemical potential 
{fi/U = 1.05) for different tunnelling rates J, (top) in the ordered case (A = 0) and (bottom) in 
the disordered case (A/t/ ~ 0.6) as predicted by the LMF (blue) and SMF (red) theory. The crosses 
(x) at the ^/j-axis represent the mean of the PD, which is the average SF order parameter -0. In the 
ordered case the PD is a delta function. With disorder all local SF parameters are zero for the MI; 
the PD in the BG phase is a superposition of a delta function at -0 = and a SF tail; the SF phase is 
characterized by a broad distribution of positive non-zero local SF parameters. 



This situation changes if disorder is introduced, since then the on-site energy is 
different on every site resulting in a variety of different values of the local SF parameter 
ijji. In the MI regime the PD is a sharp delta function at 0; = still and becomes a broad 
distribution in the BG and SF phase. In the BG phase sites with zero local SF parameter 
(corresponding to MI sites with Si = 0) coexist with sites, which have non-zero local SF 
parameter (corresponding to SF sites with 5"^ = 1) and where called SF islands before. 
In the SF regime the PD is a broad distribution representing the variety of positive 
values for the local SF parameter. Due to this characteristic behaviour the PD can 
be written as a superposition of a delta function at ^ = and a broad distribution 
representing the values ip > 0: 

P{i:) = a6{ij) + PsFW. (20) 

We denote this distribution Psf W, since it represents sites, which we refereed to as SF 
sites {Si = 1) before: 
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Starting with the same Hamihonian (11) as LMF, SMF theory introduces the 
additional approximation (15) yielding in a self-consistent equation for the PD itself, 
which is given by equation (12). The resulting distributions are shown in Figure 6 
in red. Whereas no deviations between LMF and SMF theory occur in the ordered 
case (A = 0), they become visible in the disordered case. For A > oth distribution 
have similar shapea, but especially for small values of ip, which are crucial for the 
determination of the phase transition, they differ significantly in the BG and SF phase. 

In order to identify the reason for these deviations we scrutinized the validity of 
the additional SMF restriction (15). We checked its validity by comparing the LMF 
results for the product of the PD of two different sites P (ipi) P {ipj) with the pair PD 
T^z {i'ii i'j) and determined their deviation Ap in Appendix A. As shown in Figure Al 
the approximation (15) is best in the MI regime and the BG for very small tunnelling 
rates. But for increasing JZ/U is becomes worse and especially at the phase boundary it 
fails. In Figure A2 where P (ipi) P (ipj) are Vz {ipi-, 4'j) ^i-re shown for different parameters 
it is visible that they disagree especially for small values of the ip. This disagreement 
is due to the presence of correlations of the local SF parameters ipi at different sites in 
the vicinity of the transition points, which are neglected in SMF theory. 

5. Conclusion 

In this paper we have introduced a new criterion to identify the different phases of the 
disordered BHM in d > 2 on the basis of the complete set of local boson occupation 
numbers {ui} of each sample and applied it to the ground states calculated using the 
LMF approximation. In the MI phase all (nj are integer, in the BG phase some of them 
are non-integer and form SF clusters in a MI background and in the SF phase at least 
one of these clusters percolates. The emergence of SF clusters, with an average lateral 
size that is expected to be of the order of the SF correlation length, have a finite density, 
which gives rise to a non-vanishing mean of the average SF parameters although the 
system is not SF. The latter happens only when these SF clusters percolate, which is the 
hallmark of the BG-SF transition. Moreover, the SF clusters have a fluctuating boson 
occupation number resulting in a small but non-vanishing compressibility. Thus their 
appearance is the indicator of the MI-BG transition, i.e. from the incompressible (k = 0) 
to the compressible (k > 0) phase. Consequently, the BG phase displays arbitrarily 
small but non-vanishing values for ip and k and all approaches to determine the LMF 
phase boundaries of the disordered system on the basis of the site and disorder averaged 
parameters, like the SF order parameter ip, the compressibility n, the SF or condensate 
fraction, overestimates the SF and MI phases substantially and are doomed to fail: the 
putative phase boundaries move systematically and substantially when increasing or 
decreasing the threshold only by a small amount. 

The resulting cluster analysis phase diagram for a fixed commensurate density n = 1 
is in excellent agreement with the prediction of QMC simulations, not only qualitatively 
in reproducing the characteristic shape of the SF region in the (A-[/)-diagram, but 
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also quantitatively within the numerical error bars. This is remarkably, since other 
LMF approaches using averaged quantities, like the mean SF order parameter or the 
compressibility as indicators predict much larger MI regions in the phase diagram or 
even fail to identify the BG transition, since the used indicator varies smoothly at the 
expected phase transition. Small deviations between QMC calculations and our LMF 
cluster analysis might be due to the fact that the local occupation numbers calculated 
by using the LMF approximation deviate in some regions of the phase diagram from 
the exact expectation values. Obviously it would be desirable to calculate the latter by 
QMC simulations and to perform the cluster analysis we propose to these data. 

The questions that immediately arises in this context is: 1) Is the MI-BG transition 
in the disordered BHM exactly where SF sites occur? And more interestingly 2) Is the 
BG-SF transition actually identical with the percolation transition of SF clusters? 

Concerning question 1): Although not proven rigorously, the MI-BG transition is 
supposed to occur, where the gap Eg/2, i.e. the energy for particle-hole excitations, of 
the pure, ordered BHM is equal to the disorder strength A [2]. It seems plausible that 
when this happens individual sites or small clusters will occur, where the addition or 
removal of a particle does not cost energy and thus the local boson occupation number 
fluctuates, i.e. (nj is non-integer. This is how we propose to identify the MI-BG 
transition. 

Concerning question 2) we argued in section 3.2 in the basis of the BHM to quantum 
rotor models that the SF stiffness will always vanish as long as SF clusters with a 
MI background do not percolate. This BG situation then is reminiscent of a c? + 1- 
dimensional, classical XY model with columnar disorder, in a state with (quasi) phase 
ordered finite clusters in a phase disordered background. Application of a phase twist at 
the system boundaries will only cost a macroscopic amount of energy when the ordered 
regions actually percolate - in the SF region. Note that this argument is based upon 
the existence of true long-range order in the SF phase of the pure, ordered BHM, thus 
we expect it to be valid for d >2. 

A complementary picture is based on the path-integral computation of the SF 
density [43], which is used in QMC simulations to identify SF order. The SF density or 
stiffness is proportional to the mean-square of the winding number of boson world lines 
in the path integral representation. When on average a finite fraction of boson world 
lines wrap around the whole system, the mean-square winding number is positive and the 
system is SF. To wrap around the whole system (with periodic boundary conditions), 
a boson world line, on its way through imaginary time, has to move along a path 
that traverses the whole system, thus attributing particle number fluctuations to the 
individual sites of this path. These sites will consequently attain non-integer expectation 
values for the boson occupation numbers (n-) (since for some time the boson was there 
and for some time not), thus in the end there must be at least one percolating SF cluster 
in the system. 

It should be noted that other quantum phase transitions of disordered systems 
are naturally percolation transitions, too: The critical point of the random transverse 
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Ising model is governed by an infinite randomness fixed point (in rf > 1 dimensions 
[30, 3, 44, 45]), which signals the percolation of strongly coupled clusters that away from 
criticality constitute the Griffiths phase. The percolation transition that we observe in 
our calculations falls into the universality class of a conventional, 2d site percolation 
- which means it does not carry the signature of the critical properties of the proper 
BG-SF transition. This is most probably a consequence of the LMF approximation that 
we use, since it does not properly account for spatial correlations - if applied to the 
exact ground state one would expect the critical exponents of the percolation transition 
to be related to the critical exponents of the BG-SF transition. 

In addition to providing an intuitive picture and a deeper understanding of the 
underlying physics of the phase transitions in the BHM the cluster analysis may serve 
as a reliable tool to locate the transitions in situations, in which the application of other 
criteria to discriminate the different phase might lead to erroneous predictions - as for 
instance in LMF theories. Finally, since experiments recently reached the regime of 
single site detection [9, 10] and are now able to observe the particle numbers at each 
site, an experimental application of the cluster analysis that we propose appears in 
reach. 
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Appendix A. Local SF order parameter correlations 

In this appendix we clieck tlie statistical independence assumption underlying SMF 
theory. Additional to the LMF approximation (4) made in the tunnelling part in 
the Hamiltonian, SMF theory assumes that the local SF parameters ipi, . . . ,ipz of the 
Z neighbours of a chosen site i are uncorrelated and identically distributed which is 
introduced by the approximation (15). On the basis of LMF calculations we want 
to test this approximation by comparing P (tpi) P (ipj) and Vz {ipi-, 'ipj), of which some 
examples are shown in Figure A2. The function P (tpi) P (tpj) is the product of the PD, 
describing the distribution of the local SF parameter as discussed in section 4.3. The 
function Vz {4'iy4'j) is the PD of pairs {ipi,ipj), where i and j are neighbouring sites. 
It represents the probability of having a specific value for the pair {ipi,ipj). Exactly as 
P (ipj) the Vz i'ipi, 'ipj) is computed for every realization of disorder and finally averaged. 
Both distributions should coincide if the assumption (15) is valid. In Figure Al the 
integral difference 



Ap = J dij^j d^lj,\P{ij,)Pi^lj,)-VziA,i^j)\ (A.l) 

of both distributions is shown in parameter space. 

In the MI region, where the P [ip) is a delta function at ?/^ = and for very small 
tunnelling rate JZ/U the deviations are small, whereas they are significant in the region 
of the phase transition and in the SF regime. For illustration both PDs are shown in 
Figure A2 along a line of ^/U = 1.0455 and at the tip of the Mott lobe, where the 
deviation Ap reaches its maximal value (corresponding to the black dots in Figure 
Al). Additional to the fact that all distributions are symmetric naturally, P (ipi) P (ipj) 
shows a rectangular symmetry, which intrinsically follows from the fact that it is a 
product of the same PD P (ip). The PD Vz i'^i'^j) containing further information of the 
occurring pairs shows systematic deviations. Whereas the values on the diagonal are 
reproduced quite well, the off-diagonal contributions are squeezed to the diagonal. This 
is especially pronounced in Figure A2 (d) and (h), which corresponds to the tip of the 
Mott lobe. These mean differences in comparison with P (ipi) P (ipj) can be observed 
for all parameters shown in Figure A2 and mainly occur in the regime of small local 
SF parameters. Figure Al illustrates that the assumption (15) made in SMF theory, is 
well fulfilled in the MI regime but becomes worse in the region of the phase transition. 
Whether this theory predicts the phase transition reliably in this regime, is therefore to 
question. 

The deviations occurring for small SF parameters in the limit of small ip as shown 
in Figure 6 have also been discussed in [24]. In this work the authors concluded that 
LMF theory overestimates the phase coherence in the BG regime. But this is also true 
for SMF theory, since it is based on the same approximation of the tunnelling term of 
the Hamiltonian. In this paper we resolved this apparent problem by interpreting the 
SF-BG phase transition as a percolation transition (c.f. section 4.1 and 4.2). 
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Figure Al: The deviation Ap between P {tpi) P {ijjj) and Vz (ipiii^j) is shown in dependence of the 
system parameter. In the MI regime and for very small tunnelling rates the deviations are small, 
whereas they grow at the phase transitions and in the SF regime. The black dots make the parameters 
used in Figure A2 along a line fi/U — 1.0455 and at the tip of the Mott lobes, where the deviation is 
maximal. 



P(^:,-V) 




Figure A2: The first column shows Vz {''Pijipj) and the second P {ijji) P {4'j) in the disorder case 
(A/[/ = 0.6). In the first row the parameters are given by JZ/U = 0.0283, ji/U = 1.0455 followed by 
JZ/U = 0.0586, fi/U = 1.0455 and JZ/U = 0.1414, fi/U = 1.0455 and JZ/U = 0.1434, jj/U = 0.4394 
in the last row corresponding the black dots in Figure Al. 



Appendix B. The characteristic shapes of the PD 

In section 4.3 we discussed three different shapes of the PD in the disordered case, which 
are depicted at the bottom of Figure 6. In the MI phase the PD is given by a delta 
function at ^/^ = 0, whereas in the BG and SF phase a broad distribution occurs, which 
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means that P {ip) can be represented by a superposition of a delta function at ip = 
and a continuous part Psf (V^) caused by SF sites with ^ > as defined in equation 
(20), and in the following denoted as SF distribution. Here, we will identify regions of 
the three different shapes of P {ip) in parameter space and discuss their connection to 
the phase transitions determined in section 4.1 and 4.2. 
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Figure Bl: Regions of tlic tlircc cliaractcristic sliapc of P('0) for fixed density (fi) = 1 (left) and 
fixed disorder A/[/ = 0.6 (right). The black line represents the MI-BG transition according to the 
perturbative result given by (17). Inside of the blue line P (■0) consists only of a delta function at 
■0 = 0. Within the region bounded by the red line on left and to the right of the red line on the right 
P(?/') only has a continuous part Psf (0). In other parts of the parameter space is is a superposition 
of the delta function at -0 = and the Psf ('0)- 



Numerically we identify two different benchmarks of the histograms representing 
P (ip) that we generate: The first one is the value of the histogram at the first bin, Pq, 
representing the potential delta peak of P{ip)- The second characteristic point is the 
value of the histogram at the second bin. Pi, which is given by Pi = P (■?/' = 5^) with 
5^ being the bin size of the histogram representing P {if)) [5^ = 0.0025 in LMF and 
5^ = 0.015 in SMF theory). 

In Figure Bl the regions determined on the basis of these benchmarks are shown 
for fixed density n = 1 on the left and fixed disorder strength A/U = 0.6 on the right. 
The LMF results are depicted with circles, SMF results with crosses. The blue curves 
enclose the regions in which P [tp) is just a delta function at ^/^ = (numerically: Pq > 
and Pi < 10~^), i.e. the system contains only MI sites and is therefore in the MI phase. 
The red curves delimit the regions, in which the delta function part of P {ip) vanishes 
(numerically Pq < 10~^), which means that all sites are SF sites. In the remaining part 
of the phase diagram P (tp) consists of a superposition of a delta function at ?/^ = 
and a continuous part Psf (^), i.e. the system has MI and SF sites. In this region the 
system can either be in the BG phase, or, if the SF sites percolate, in the SF phase. 
Therefore only the MI-BG phase boundary can be extracted from the characteristics 
of the probability distribution of the local order parameter, P (ip), but not the BG-SF 
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boundary. 

LMF and SMF theory coincide very well at the blue line, which describes the 
occurrence of the SF distribution. Deviations are visible at the red line, where the 
delta function at ^ = disappears and the PD is purely given by Psf W- While these 
discrepancies are rather small for small disorder strengths and U/J > 23, they enlarge 
with increasing disorder. In Appendix A we tested the validity the SMF approximation 
(15) by comparing the LMF results for the product of the PD of two different sites 
P (ipi) P (ipj) with the pair PD Vz (ipi, ipj)- In Figure Al its deviation Ap is shown in 
dependence of fi/U and JZ/U for fixed disorder A/U = 0.6. In comparison with the 
diagram on the right in Figure Bl it is obvious that in the region of the blue line the SMF 
assumption (15) for P (tp) is valid. However, at the red curve deviations between LMF 
and SMF theory occur, since close to the BG-SF phase transition the approximation 
(15) is expected to be invalid. 
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